Statistics Project Visualizations

library(readxl)
library(dplyr)
library(ggplot2)
library(knitr)
library(tidyr)
library(sf)
library(plotly)

data <- read_excel("data.xlsx", sheet = "LFS_Year_2023")

regions <- c(
  "Yerevan",
  "Aragatsotn",
  "Ararat",
  "Armavir",
  "Gegharkunik",
  "Lori",
  "Kotayq",
  "Shirak",
  "Syuniq",
  "Vayoc Dzor",
  "Tavush"
)

people_by_regions <- data %>%
  group_by(A3) %>%
  summarise(
    n = n()
  ) %>%
  mutate(A3 = factor(A3, levels = 1:11, labels = regions)) %>%
  rename(
    region = A3,
    population = n
  )

people_by_regions

ggplot(people_by_regions, aes(x = region, y = population)) +
  geom_col(fill = "plum3") +
  theme_minimal() +
  labs(
    title = "Population by Regions",
    x = "Regions",
    y = "Population"
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

urban_rural <- data %>%
  group_by(A5) %>%
  summarise(n = n()) %>%
  mutate(A5 = factor(A5, levels = c(1, 2), labels = c("Urban", "Rural")))

urban_rural

gender_ratio <- data %>%
  group_by(B3) %>%
  summarise(n = n()) %>%
  mutate(B3 = factor(B3, levels = c(1, 2), labels = c("Male", "Female")))

gender_ratio

age_groups <- c(
  "0 - 4",
  "5 - 9",
  "10 - 14",
  "15 - 19",
  "20 - 24",
  "25 - 29",
  "30 - 34",
  "35 - 39",
  "40 - 44",
  "45 - 49",
  "50 - 54",
  "55 - 59",
  "60 - 64",
  "65 - 69",
  "70 - 74",
  "75+"
)

people_by_age <- data %>%
  group_by(Age_16groups) %>%
  summarise(n = n()) %>%
  mutate(Age_16groups = factor(Age_16groups, levels = 1:16, labels = age_groups)) %>%
  rename(
    age_group = Age_16groups,
    population = n
  )

people_by_age

ggplot(people_by_age, aes(x = age_groups, y = population)) +
  geom_col(fill = "plum3") +
  theme_minimal() +
  labs(
    title = "Population by Age Group",
    x = "Age Group",
    y = "Population"
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

educational_levels <- c("Illiterate", "No primary", "Primary", "Basic", "Secondary / high", 
                      "Vocational", "Secondary specialized", "Bachelor's degree", 
                      "Master's degree", "Certified specialist", 
                      "Post-graduate (Ph.D, doctorate)")

people_by_education <- data %>%
  group_by(B7) %>%
  summarise(n = n()) %>%
  mutate(B7 = factor(B7, levels = 1:11, labels = educational_levels)) %>%
  rename(
    educational_level = B7,
    population = n
  )

people_by_education

ggplot(people_by_education, aes(x = educational_level, y = population)) +
  geom_col(fill = "plum3") +
  theme_minimal() +
  labs(
    title = "Population by Educational Level",
    x = "Educational Level",
    y = "Population"
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

nace_sectors <- c(
  "Agriculture, Forestry and Fishing",
  "Mining and Quarrying",
  "Manufacturing",
  "Electricity, Gas, Steam and Air Conditioning Supply",
  "Water Supply; Sewerage, Waste Management and Remediation Activities",
  "Construction",
  "Wholesale and Retail Trade; Repair of Motor Vehicles and Motorcycles",
  "Transportation and Storage",
  "Accommodation and Food Service Activities",
  "Publishing, Broadcasting, and Content Production and Distribution Activities",
  "Telecommunication, Computer Programming, Consulting, Computing Infrastructure, and other Information Service Activities",
  "Financial and Insurance Activities",
  "Real Estate Activities",
  "Professional, Scientific and Technical Activities",
  "Administrative and Support Service Activities",
  "Public Administration and Defence; Compulsory Social Security",
  "Education",
  "Human Health and Social Work Activities",
  "Arts, Sports and Recreation",
  "Other Service Activities",
  "Activities of Households as Employers; Undifferentiated Goods and Services Producing Activities of Households for Own Use"
)

employees_by_industry <- data %>%
  group_by(E4_21groups_NACE_rev_2.2) %>%
  summarise(n = n()) %>%
  mutate(E4_21groups_NACE_rev_2.2 = factor(E4_21groups_NACE_rev_2.2, levels = 1:21, labels = nace_sectors)) %>%
  rename(
    industry = E4_21groups_NACE_rev_2.2,
    employee_number = n
  )

employees_by_industry

exact_salary_data_points <- data %>%
  select(E14_1) %>%
  na.omit() %>%
  summarise(n = n())

exact_salary_data_points

exact_salary_desc <- data %>%
  select(E14_1) %>%
  na.omit() %>%
  summarise(mean = mean(E14_1), max = max(E14_1), min = min(E14_1), median = median(E14_1))

exact_salary_desc

mean_salary_by_gender <- data %>%
  select(E14_1, B3) %>%
  mutate(B3 = factor(B3, levels = 1:2, labels = c("Male", "Female"))) %>%
  na.omit() %>%
  group_by(B3) %>%
  summarise(mean = mean(E14_1), median = median(E14_1), min = min(E14_1), max = max(E14_1))

mean_salary_by_gender

salary_data <- data %>%
  select(salary = E14_1, gender_code = B3) %>%
  mutate(
    gender = factor(gender_code, levels = c(1, 2), labels = c("Male", "Female"))
  ) %>%
  filter(!is.na(salary), !is.na(gender))

ggplot(salary_data, aes(x = gender, y = salary, fill = gender)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 500000)) +
  theme_minimal() +
  scale_fill_manual(values = c("Male" = "#56B4E9", "Female" = "#E78AC3")) +
  labs(
    title = "Salary Distribution by Gender",
    x = "Gender",
    y = "Salary"
  ) +
  theme(legend.position = "none")

income_ranges <- c(
  "Up to 55 000 AMD",
  "55 000 AMD",
  "55 001 - 110 000 AMD",
  "110 001 - 220 000 AMD",
  "220 001 - 440 000 AMD",
  "440 001 - 600 000 AMD",
  "600 001 - 700 000 AMD",
  "Refused to answer",
  "Do not know / difficult to answer"
)

categorical_salary_data_points <- data %>%
  select(E15) %>%
  na.omit() %>%
  summarise(number_of_data_points = n())

categorical_salary_data_points

categorical_salary_frequency <- data %>%
  group_by(E15) %>%
  summarise(n = n()) %>%
  na.omit() %>%
  mutate(E15 = factor(E15, levels = 1:9, labels = income_ranges)) %>%
  rename(
    income_range = E15,
    frequency = n
  )

categorical_salary_frequency

ggplot(categorical_salary_frequency, aes(x = income_range, y = frequency)) +
  geom_col(fill = "plum3") +
  theme_minimal() +
  labs(
    title = "Salary Histogram",
    x = "Income Ranges",
    y = "Frequency"
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Data Cleaning

checking to see how many rows have more than 50% missing data (from the columns we are going to work on)

cols_to_work <- data[, c('A3', 'A5', 'A6_Month', 'A6_Year', 'B3', 'Age_16groups', 'B7', 'B8', 
                        'B8_Year', 'B9', 'B10', 'B11', 'C1', 'C9', 'C4', 'C3', 'D1', 'D2', 'D3', 'E2_9group', 
                        'E4_21groups_NACE_rev_2.2', 'E5', 'E6', 'E7', 'E9', 'E10', 'E13', 'E14_1', 'E14_2', 'E15', 'E16', 
                        'E17', 'E19', 'E23', 'E24', 'F1', 'F2_9group', 'F3_21groups_NACE_rev_2.2', 'F4', 'F5', 
                        'F6', 'F8_1', 'F8_2', 'G2_1', 'G2_2', 'G2_3Total', 'G3', 'J1_4group', 'J3_1', 'J3_2', 'J3_LastJob')]
missing_per_row <- rowSums(is.na(cols_to_work)) 
total_cols <- length(cols_to_work)
exceeds <- which(missing_per_row > 0.5 * total_cols)
good_rows <-length(exceeds)
print(good_rows)

[1] 16897

check how many empty rows there are

sum(apply(data, 1, function(row) all(is.na(row))))

[1] 0

check how many rows are missing values of important columns

cols_to_check <- c('A3', 'A5', 'A6_Month', 'A6_Year', 'B3', 'Age_16groups', 'B7', 'B8', 
                        'B8_Year', 'B9', 'B10', 'B11', 'C1', 'C9', 'C4', 'C3', 'D1', 'D2', 'D3', 'E2_9group', 
                        'E4_21groups_NACE_rev_2.2', 'E5', 'E6', 'E7', 'E9', 'E10', 'E13', 'E14_1', 'E14_2', 'E15', 'E16', 
                        'E17', 'E19', 'E23', 'E24', 'F1', 'F2_9group', 'F3_21groups_NACE_rev_2.2', 'F4', 'F5', 
                        'F6', 'F8_1', 'F8_2', 'G2_1', 'G2_2', 'G2_3Total', 'G3', 'J1_4group', 'J3_1', 'J3_2', 'J3_LastJob') 
sum(rowSums(is.na(data[cols_to_check])) > 0)

[1] 27007

sum(!is.na(data$E14_2))

[1] 6532

get_mode <- function(x) {
  uniqx <- unique(x)
  uniqx[which.max(tabulate(match(x, uniqx)))]
}

data %>%
  filter(!is.na(E14_2)) %>%
  select(E14_2) %>%
  summarise(mode = get_mode(E14_2))

data %>%
  filter(!is.na(E14_2)) %>%
  count(E14_2, sort = TRUE)

Descriptive Data Analysis and Visualizations

cleaned_data <- read_excel("cleaned_data.xlsx")

Barplot: Counts by Age Range and Gender

Decoding data

decoded_data <- cleaned_data %>%
  mutate(
    B3 = case_when(
      B3 == 1 ~ "Male",
      B3 == 2 ~ "Female",
      TRUE ~ NA_character_
    )
  ) %>%
  mutate(
    Age_16groups = case_when(
      Age_16groups == 1 ~ "0-4",
      Age_16groups == 2 ~ "5-9",
      Age_16groups == 3 ~ "10-14",
      Age_16groups == 4 ~ "15-19",
      Age_16groups == 5 ~ "20-24",
      Age_16groups == 6 ~ "25-29",
      Age_16groups == 7 ~ "30-34",
      Age_16groups == 8 ~ "35-39",
      Age_16groups == 9 ~ "40-44",
      Age_16groups == 10 ~ "45-49",
      Age_16groups == 11 ~ "50-54",
      Age_16groups == 12 ~ "55-59",
      Age_16groups == 13 ~ "60-64",
      Age_16groups == 14 ~ "65-69",
      Age_16groups == 15 ~ "70-74",
      Age_16groups == 16 ~ "75+",
      TRUE ~ NA_character_
    )
  ) %>%
  mutate(
    E15 = case_when(
      E15 == 1 ~ "Up to 55 000 AMD",
      E15 == 2 ~ "55 000 AMD",
      E15 == 3 ~ "55 001 -  110 000 AMD",
      E15 == 4 ~ "110 001 - 220 000 AMD",
      E15 == 5 ~ "220 001 - 440 000 AMD",
      E15 == 6 ~ "440 001 - 600 000 AMD",
      E15 == 7 ~ "600 001 - 700 000 AMD",
      E15 == 8 ~ "Refused to answer",
      E15 == 9 ~ "Do not know / difficult to answer",
      TRUE ~ NA_character_
    )
  ) %>%
  mutate(
    A3 = case_when(
      A3 == 1 ~ "Yerevan",
      A3 == 2 ~ "Aragatsotn",
      A3 == 3 ~ "Ararat",
      A3 == 4 ~ "Armavir",
      A3 == 5 ~ "Gegharkunik",
      A3 == 6 ~ "Lori",
      A3 == 7 ~ "Kotayq",
      A3 == 8 ~ "Shirak",
      A3 == 9 ~ "Syuniq",
      A3 == 10 ~ "Vayoc Dzor",
      A3 == 11 ~ "Tavush",
      TRUE ~ NA_character_
    )
  ) %>%
  mutate(
    A5 = case_when(
      A5 == 1 ~ "Urban",
      A5 == 2 ~ "Rural",
      TRUE ~ NA_character_
    )
  )

ggplot(decoded_data, aes(x = factor(Age_16groups), fill = factor(B3))) +
  geom_bar(position = "stack") + 
  labs(
    title = "Age Distribution by Gender",
    x = "Age Range",
    y = "Count",
    fill = "Gender"
  ) +
  scale_fill_manual(values = c("Male" = "#56B4E9", "Female" = "#E78AC3")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

Also, let’s calculate this numerically.

gender_counts <- decoded_data %>%
  group_by(Age_16groups, B3) %>%
  summarise(Count = n(), .groups = "drop")

gender_counts_wide <- gender_counts %>%
  pivot_wider(names_from = B3, values_from = Count, values_fill = 0) %>%
  rename("Age Range" = Age_16groups)

# display
knitr::kable(gender_counts_wide, caption = "Gender Count by Age Group")

Gender Count by Age Group

Age Range	Female	Male
15-19	29	56
20-24	177	322
25-29	273	576
30-34	418	749
35-39	552	762
40-44	557	626
45-49	492	540
50-54	444	481
55-59	456	491
60-64	535	585
65-69	282	341
70-74	105	159
75+	44	58

Findings: In this chart, we have considered the records where the implicit salaries are mentioned. Both the barchart and the tables highlight that the number of males taking the survey is higher than the number of women, independent of their age range.

Bar Plot: Mean Salary by Age Range and Gender

age_gender_stats <- decoded_data %>%
  filter(!is.na(E14_1), !is.na(Age_16groups), !is.na(B3)) %>%
  group_by(Age_16groups, B3) %>%
  summarize(
    mean_salary = mean(E14_1),
    se = sd(E14_1) / sqrt(n()),
    lower_ci = mean_salary - 1.96 * se,
    upper_ci = mean_salary + 1.96 * se,
    n = n(),
    .groups = "drop"
  ) %>%
  filter(n >= 10) %>%
  mutate(salary_label = paste0(format(round(mean_salary), big.mark = ",")))

ggplot(age_gender_stats, 
       aes(x = Age_16groups, y = mean_salary, 
           color = B3, group = B3)) +
  geom_line(linewidth = 1.2) +
  geom_pointrange(
    aes(ymin = lower_ci, ymax = upper_ci),
    position = position_dodge(width = 0.2),
    size = 0.7
  ) +
  geom_text(aes(y = upper_ci, label = salary_label),
            position = position_dodge(width = 0.5),
            vjust = -0.5, size = 3, color = "black") +
  labs(
    title = "Mean Salary by Age Group and Gender with 95% Confidence Intervals",
    x = "Age Category", 
    y = "Mean Salary (AMD)",
    color = "Gender"
  ) +
  scale_y_continuous(labels = scales::comma) +
  theme_minimal()

Findings: In this chart, we have considered the records where the implicit salaries are mentioned. We can see that the salary of males in the sample is considerably higher and the condifence interval is bigger. However, the length of confidence intervals for the mean salary of males are longer, which means that we are not as much certain about the mean as for the females. The lower bounds of the confidence intervals of males’ mean salaries are closers to the females’ mean salaries.

Mean Salary with 95% Confidence Intervals

decoded_data <- decoded_data%>%
  mutate(
    B8 = case_when(
      B8 == 1 ~ "Up tp 2000",
      B8 == 2 ~ "2001-2005",
      B8 == 3 ~ "2006-2010",
      B8 == 4 ~ "Since 2012"
    )
  )

Based on the salary they have provided.

correct_order <- c("Up tp 2000", "2001-2005", "2006-2010", "Since 2012")

plot_data <- decoded_data %>%
  mutate(
    graduation_year_ordered = factor(
      B8,
      levels = correct_order,
      ordered = TRUE
    )
  ) %>%
  filter(!is.na(E14_1)) %>%
  group_by(graduation_year_ordered, B3) %>%  # Group by the ordered factor
  summarize(
    mean_salary = mean(E14_1),
    se = sd(E14_1) / sqrt(n()),
    lower_ci = mean_salary - 1.96 * se,
    upper_ci = mean_salary + 1.96 * se,
    .groups = "drop"  # Explicitly drop grouping
  ) %>%
  mutate(salary_label = paste0(format(round(mean_salary), big.mark = ","), " AMD"))

ggplot(plot_data, aes(x = graduation_year_ordered, y = mean_salary, 
                     color = B3, group = B3)) +
  geom_pointrange(aes(ymin = lower_ci, ymax = upper_ci),
                 position = position_dodge(width = 0.5)) +
  geom_line(position = position_dodge(width = 0.5)) +
  geom_text(aes(y = upper_ci, label = salary_label),
            position = position_dodge(width = 0.5),
            vjust = -0.5, size = 3, color = "black") +
  labs(title = "Mean Salary with 95% Confidence Intervals",
       x = "Graduation Year", y = "Mean Salary (AMD)",
       color = "Gender") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_segment()`).
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_text()`).

Let’s also see the count of men and women based on the graduation year.

graduation_summary <- decoded_data %>%
  filter(!is.na(B3), !is.na(B8)) %>%
  count(Gender = B3, Graduation_Year = B8, name = "Count") %>%

  group_by(Graduation_Year) %>%
  mutate(Percentage = Count / sum(Count) * 100) %>%
  ungroup() %>%
  
  mutate(Percentage_Label = paste0(round(Percentage, 1), "%")) %>%
  
  select(Gender, Graduation_Year, Count, Percentage, Percentage_Label) %>%
  arrange(Graduation_Year, Gender)

print(graduation_summary)

Findings: In this chart, we have considered the records where the implicit salaries are mentioned. We see that the number of females is bigger here. However, the mean salary of the males is still higher, with bigger confidence intervals as in the previous chart.

Can We Find the Mean Salary by Work Location and Gender

decoded_data <- decoded_data%>%
mutate(
    J3_1 = case_when(
      J3_1 == 1 ~ "Armenia",
      J3_1 == 2 ~ "Abroad",
      J3_1 == 3 ~ "No Paid Job"
      
    )
  ) %>%
mutate(
    J3_2 = case_when(
      J3_2 == 1 ~ "Armenia",
      J3_2 == 2 ~ "Abroad",
      J3_2 == 3 ~ "No Paid Job"
      
    )
  )

unique(cleaned_data['J3_1'])

unique(decoded_data['J3_2'])

It comes out that we cannot work with columns J3_1 and J3_2, which specify if the survey taker works in Armenia/abroad/both/no paid job. The reason is that when cleaning the data, we kept only the records where salary range or exact salary is mentioned.

decoded_data <- decoded_data%>%
  mutate(
    E10 = case_when(
      E10 == 1 ~ "Public",
      E10 == 2 ~ "NGO",
      E10 == 3 ~ "Private",
      E10 == 4 ~ "Self-emplyed"
    )
  ) %>%
  mutate(
    F6 = case_when(
      F6 == 1 ~ "Public",
      F6 == 2 ~ "NGO",
      F6 == 3 ~ "Private",
      F6 == 4 ~ "Self-employed"
    )
  )

Let’s see how many records we have with filled primary job and secondary ones.

primary_job_count <- decoded_data %>%
  filter(!is.na(E10)) %>%  
  count(B3, name = "count")  
primary_job_count

secondary_job_count <- decoded_data %>%
  filter(!is.na(F6)) %>%  
  count(B3, name = "count")

secondary_job_count

plot_data <- decoded_data %>%
  select(gender = B3, salary = E14_1, ownership = E10) %>%  
  filter(!is.na(salary), !is.na(ownership)) %>%              
  mutate(ownership = factor(ownership)) 

summary_stats <- plot_data %>%
  group_by(gender, ownership) %>%
  summarize(
    mean_salary = mean(salary),
    se = sd(salary) / sqrt(n()),
    lower_ci = mean_salary - 1.96 * se,
    upper_ci = mean_salary + 1.96 * se,
    n = n(),
    .groups = "drop"
  ) %>%
  filter(n >= 10) 

ggplot(summary_stats, 
       aes(x = ownership, y = mean_salary, fill = gender)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_errorbar(
    aes(ymin = lower_ci, ymax = upper_ci),
    position = position_dodge(width = 0.8),
    width = 0.2,
    color = "black"
  ) +
  geom_text(
    aes(label = paste0(round(mean_salary/1000, 1), "k \n(n=", n, ")"),
        y = upper_ci),
    position = position_dodge(width = 0.8),
    vjust = -0.5,
    size = 3
  ) +
  labs(
    title = "Mean Salary (AMD) by Company Ownership (Primary Job)",
    subtitle = "Error bars show 95% confidence intervals",
    x = "Ownership Type",
    y = "Mean Salary (AMD)",
    fill = "Gender"
  ) +
  scale_y_continuous(labels = scales::comma) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Findings: In this chart, we have considered the records where the implicit salaries are mentioned. Concentrating on the chart above, we can see that only females are recorded as self-employed (keeping in mind that the data is filtered to keep only records with mentioned salary). The analysis also shows that males earn more in NGOs, Private and Public enterprises, but the highest difference is in Private companies. Although there are more females who answered these questions then males, the difference still exists.

mean_salary_by_region <- decoded_data %>%
  select(gender = B3, region = A3, salary = E14_1) %>%
  filter(!is.na(region), !is.na(gender), !is.na(salary)) %>%
  group_by(region, gender) %>%
  summarise(mean_salary = mean(salary, na.rm = TRUE),
            median_salary = median(salary, na.rm = TRUE), 
            min_salary = min(salary, na.rm = TRUE), 
            max_salary = max(salary, na.rm = TRUE))

`summarise()` has grouped output by 'region'. You can override using the `.groups` argument.

  
mean_salary_by_region

difference_in_mean_salary_by_region <- mean_salary_by_region %>%
  pivot_wider(names_from = gender, values_from = mean_salary) %>%
  mutate(wage_gap = Male - Female) %>%
  select(region, wage_gap)

difference_in_mean_salary_by_region

armenia_map <- st_read("armenia.json")

Reading layer `armenia' from data source `D:\Statistics\Final Project\armenia.json' using driver `GeoJSON'
Simple feature collection with 11 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 43.43629 ymin: 38.8637 xmax: 46.60261 ymax: 41.29045
Geodetic CRS:  WGS 84

region_name_map <- tibble::tibble(
  map_name = c("Tavush", "Lori", "Shirak", "Gegharkunik", "Vayots Dzor", 
               "Syunik", "Ararat", "Aragatsotn", "Armavir", "Kotayk", "Erevan"),
  region = c("Tavush", "Lori", "Shirak", "Gegharkunik", "Vayoc Dzor",
             "Syuniq", "Ararat", "Aragatsotn", "Armavir", "Kotayq", "Yerevan")  # match your data
)

armenia_map <- armenia_map %>%
  left_join(region_name_map, by = c("name" = "map_name"))

wage_gap_map <- armenia_map %>%
  left_join(difference_in_mean_salary_by_region, by = "region")

ggplot(wage_gap_map) +
  geom_sf(aes(fill = wage_gap), color = "white") +
  scale_fill_gradientn(
    colors = c("navy", "skyblue", "yellow", "orange", "red"),
    na.value = "gray90",
    name = "Wage Gap (AMD)"
  ) +
  labs(
    title = "Gender Wage Gap by Region (Armenia)",
    subtitle = "Difference in mean salary: Male - Female",
    caption = "Data: Your Survey"
  ) +
  theme_minimal()

p <- ggplot(wage_gap_map) +
  geom_sf(aes(fill = wage_gap, text = paste("Region:", name, "<br>Wage Gap:", wage_gap))) +
  scale_fill_viridis_c() +
  theme_minimal()

Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,  :
  Ignoring unknown aesthetics: text

ggplotly(p, tooltip = "text")

Findings: As shown in the map, there is a wage gap in all regions of Armenia, which fluctuates roughly between 50000 and 25000 AMD, with the smallest gap present in Ararat and the biggest in Armavir.

Urban vs Rural mean salaries

urban_rural_salary <- decoded_data %>%
  select(settlement = A5, gender = B3, salary = E14_1) %>%
  filter(!is.na(settlement), !is.na(gender), !is.na(salary)) %>%
  group_by(settlement, gender) %>%
  summarise(mean_salary = mean(salary), min_salary = min(salary), max_salary = max(salary))

`summarise()` has grouped output by 'settlement'. You can override using the `.groups` argument.

urban_rural_salary

ggplot(decoded_data, aes(x = A5, y = E14_1, fill = A5)) +
  geom_boxplot(outlier.color = "red", outlier.alpha = 0.3) +
    coord_cartesian(ylim = c(0, 500000)) +
  labs(
    title = "Income Distribution: Urban vs Rural",
    x = "Settlement Type",
    y = "Salary (AMD)"
  ) +
  theme_minimal() +
  scale_fill_manual(values = c("Urban" = "#66c2a5", "Rural" = "#fc8d62"))

Warning: Removed 3575 rows containing non-finite outside the scale range (`stat_boxplot()`).

Findings: The mean salary in urban settlements is higher than in rural settlements. Moreover, the range of salaries is much bigger than that of urban salaries.